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ABSTRACT 

We present a new method for determining time delays among the light curves of various 
images in a gravitational lens. The method is based on constructing a simple model for the 
source variation and forming a y 2 measure of the agreement of this same variation with all 
of the lightcurves. While inspired by Press et al. (1992a, b) our approach is different since 
we do not assume a Gaussian process for the source variation. Our method has a number 
of desirable properties: first, it yields an approximate reconstruction of the source variation 
and of other parameters such as relative time delays; second, it easily incorporates different 
assumptions about the relations among the light curves and about the data measurement 
errors; finally, it can be applied to more than two light curves by addition of y 2 . We apply 
this method to the light curves of the quadruple gravitational lens PG1115+080 measured by 
Schechter et al. (1997). Unlike Schechter et al. we include correlated measurement errors in 
the analysis, as well as the possibility that microlensing may cause different images to vary 
by different factors in flux. We find a value of 25.0^31 days (95% confidence) for the delay 
between components B and C (close to the 24 day value of Schechter et al., and so leading to 
a similar value of the Hubble constant for a given lens model). However, the ratio t ac/ tBA 
of the two intermediate delays is poorly determined at 1.131^7 (68% confidence), close to the 
value predicted by lens models (~ 1.4) unlike the Schechter et al. value (~ 0.7). The variation 
ratios of C with respect to A and of A with respect to B are both different from 1, 1.39l'2Q 
and .79lt'i2 (95% confidence), respectively. This is an indication of a microlensing gradient, 
and this type of microlensing may allow us to conclude that the size of the quasar optical 
emission region is about 1000 AU. 


Subject headings: cosmology — distance scale — gravitational lensing 


Introduction 


Even before the discovery of the first gravitational lens 0957+561 ( |Walsh, Carswell &: Weymann 19~79 ) 
it was recognized that measurements of the time delay between images can be used to determine the 
Hubble constant ( Refsdal 1964 1966). Despite a history of systematic difficulties, recent measurements 
combined with an analysis of lens models (Grogin &; Narayan 1996a, 1996b|) indicate that a robust 
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measurement of the Hubble constant (Hq) with an accuracy comparable to that of more conventional 
techniques may be within reach. The measurements include a precise optical time delay (Kundic et al. 
1997) which is consistent with the latest results from radio monitoring (Haarsma et al. 1997), and a 
measurement of the velocity dispersion of the lens galaxy (Falco et al. 1997). 

Since each lens can potentially yield an independent, single step Ho measurement, obtaining values in a 
number of lenses with a variety of morphologies and constraints on models may help eliminate systematic 
errors. Each lens requires a well-constrained lens model, lens and source redshifts, and a measurement of 
the time delays among the images. Observational constraints on the lens model may include precise image 
positions and flux ratios as well as data on the objects responsible for the lensing. Resolved structure 
in images as revealed by radio interferometry provides many constraints, since it is essentially the same 
as observing multiple sources with the same lens. Whether lensing involves a galaxy, group, or cluster, 
the position, velocity dispersion, and other observational probes of the mass distribution of the lensing 
objects all yield constraints on lens models. Measuring time delays requires observing time variations in 
the image intensities with sufficient accuracy and time resolution. As the number of lenses being carefully 
monitored increases, more time delays are being determined, such as the preliminary measurement in PKS 
1830-211 ( |Van Ommen et al. 1995| ). Flux measurements are not the only possibility, as shown by the 
promising measurements of variations in polarization fraction in the images of B 0218+357 (Corbett et 
al. |l996| ). 


The quadruply imaged quasar PG1115+080 was the second lens to be discovered (Weymann et al 


198(0). It is radio quiet, but optical Hubble Space Telescope images ( Kristian et al. 1993|) were recently 
analyzed by Schechter et al. (1997, hereafter SCH) and by Keeton & Kochanek (1997). They found that 
lens models which include the effect of the lens galaxy and that of the nearby group of galaxies discovered 
by Young et al. (1981) can fit the image positions well. They however still found great freedom in the 
Ho values predicted by these lens models for a given time delay. 

In four-image configurations where the images lie at roughly the same distance from the lens, there is 
a well known degeneracy between the radial profile of the lens mass and the inferred Hq (Kochanek 1991; 


Wambsganss fe Paczyriski 1994 ). In this situation the relative image positions do determine the total 
enclosed mass within the ring of images, but they are not very sensitive to the radial profile of the mass. 
Changing the radial profile affects the convergence at the images, and this changes the inferred Ho value 
in a similar way to the effect of a constant mass sheet ( Falco et al. 1985 ; Narayan 1991 ). In PG1115+080 
future observations of the lensing galaxy light profile and, ultimately, a direct measurement of its central 
velocity dispersion may constrain or eliminate this degeneracy. 


Recently SCH measured light curves for the Al, A2, B and C images of PG1115+080, and used them 
to determine multiple delays. The bright Al and A2 images are close together and should have a very 
small time delay (~ a few hours), so they were combined into a single A=A1+A2 curve. SCH used the 
Press et al. (1992a, b) method and found that C leads A and A leads B, with tBC = 23.7 ± 3.4 days and 
tAC = 9-4 ± 3.4 days which yields a ratio tabc = tAc/^BA of 0.7 ± 0.3. It is useful to express the two 
independent quantities as tsc and tabc, since tsc can be taken to fix the Hubble constant for a given 
lens model while i'abc is independent of overall distance, and can be compared directly with the ratio 
predicted by lens models. The models mentioned above are consistent in predicting tabc = 1-4 to within 
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about 0.1, and SCH noted the 2a discrepancy with their fitted value. In their analysis SCH assumed that 
the measurement errors in the light curves are uncorrelated, and also that the fractional flux variations 
are the same for each component. 

In this paper we present a more detailed analysis of the light curves in PG1115+080. We first present 
a new method based on y 2 minimization, which has many of the advantages of Press et al. (1992a, b), 
but is simpler and allows for a more conservative assessment of errors in the reconstructed parameters. 
We then apply this method to PG1115+080, and include correlated measurement errors in the analysis. 
In addition to relative time delays, we also allow for different factors of variation in flux, which may arise 
from microlensing of the images. In §2 we present our y 2 method, and discuss its distinct features and 
free parameters. We then discuss the physical meaning of the various parameters that the method can 
account for and attempt to extract from data. In §3 we apply our method to fitting the PG1115+080 
light curves, singly, in pairs, and all together, and discuss the results and implications. Finally in §4 we 
summarize our results and point out some of the significant returns possible from further monitoring. 


2. Method and physical parameters 

Suppose we have a light curve of N points, (ij, Vi) i = 1,.., N, where t, is time (in days) and Vi is the 
intensity (in magnitudes) at time t*, and we construct a model for v(t) with M adjustable parameters a,j, 
j = 1,..., M. Then if the Vi measurements are independent, with normally distributed errors of standard 
deviation a t for each i. the y 2 measure of goodness-of-fit is 


N 


x 2 = E 


i= 1 L 


v(ti, 0 ,1 . . . CIm) 


£T; 


( 1 ) 


When the measurement errors are correlated, with a covariance matrix Q, then 

N 

x 2 = E W ij \- Vi ~ V ( ti ’ °1 • • ■ ®m)] [Vj - v(tj\a 1 ... a M )\ , (2) 

i,j= 1 

where the weight matrix W is the matrix inverse of Q. 

We wish our model for v(t) to be an arbitrary continuous curve. Although quasar variations may 
be stochastic, the variations typically decrease in amplitude at smaller and smaller time scales, and a 
continuous curve can describe the large time-scale variations well. Of course, we cannot take an arbitrary 
curve, since we must describe the curve with a small number of parameters. We therefore take L points 
(tfc, Vk) with the just covering the range of the data and equally spaced. The vk are then the parameters, 
and the model for v(t) is a curve going through all the ( tk,Vk ) with interpolation used for other values of 
t. Note that we interpolate between the Vk in the model, only as a way of parameterizing a general curve 
v(t) in terms of the L parameters; we do not interpolate between the noisy data points directly. When 
we fit multiple light curves, a single such model is used for all of them, or equivalently the light curves 
are all combined into a single light curve which is then compared to the model. Additional parameters 
can be added to account for relative shifts of an entire light curve relative to the others, before comparing 
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to the model. A relative magnification corresponds (in magnitudes) to adding a constant to all the 
measured V{ of an image. A relative time delay corresponds to adding a constant to all the measured f, 
of an image. For given time delays, all the other parameters are linear (i.e. the x 2 is a quadratic form 
in these parameters) if we use only interpolation which is linear in the Vk (but arbitrarily non-linear in 
the tfc). Minimizing the x 2 then corresponds to solving a linear system, and we solve it using Cholesky 
decomposition (e.g. Press et al. 1992cj , §2.9). 

Like Press et al. (1992a, b), our method compares a model for the source variation to the light curves 
with a x 2 measure, and thus both methods assume that the statistics of the measurement errors are 
Gaussian. However, the two methods construct the model of the source variation very differently. In 
Press et al. the source variation is assumed to be a stationary Gaussian process, with a correlation 
function which is extracted from the light curves themselves. In that method there is some freedom to 
vary the assumed correlation function in order to check for robustness. In our method we assume no 
particular form for the correlations among source fluxes at different source times. By taking a continuous 
curve we are effectively smoothing the actual source variation over some small time scale, but the curve 
is otherwise arbitrary. Although the two methods use different assumptions, they may not necessarily 
yield very different results, since the intrinsic variability of quasars appears to decrease at small time 
scales. Qualitatively, then, both methods are dominated by requiring smoothness for nearby points when 
all the light curves are combined. We also note that the correlations which we are including are in the 
observational errors of different images observed at a single time. Such correlations were not accounted 
for by SCH, but they in principle can also be included in the Press et al. method. 


Our model is very flexible. There is freedom to choose both the number of interpolation points L and 
the interpolation method. We use several methods of interpolation. Given some time t between (ffc,Vfc) 
and (ffc + i,Ufc + i), in linear interpolation we take for v(t) the linearly interpolated value between v & and 
Vk+i- In cubic interpolation we interpolate using the cubic polynomial through the four points at k — 1, 
k, k + 1 and k + 2 (with linear interpolation used for k = 1 or k = L — 1). In cubic spline interpolation 
(e.g. Press et al. 1992c] , §3.3), we draw a cubic polynomial between every pair of points k and k + 1 
so that the entire curve is continuous in the second derivative. We find in §3 below that changing the 
interpolation method has a rather minor effect. Note that with any of these interpolation methods there 
is some inter-dependence among all of the Vk, unless there is some ko with no data points between tk 0 
and ffc 0+ 1 - In this case the parameters Vk for k > ko would decouple from those with k < ko, at least with 
linear interpolation. There are 37 points in each light curve in the PG1115+080 data, and the final data 
point follows a gap of 38 days. Hereafter we drop this last point, and then we are never in a situation 
where such a ko exists (The above argument suggests that the last point has only a small effect on the 
results with our method, and this is indeed the case). 


The number of points in each light curve is thus 36, and we use an L of 10, 20 and 30. Below we 
find that the differences in the derived parameters using the various L are fairly large, so this freedom 
provides a more conservative and robust estimate of parameter uncertainties. Using small L corresponds 
essentially to a larger-scale smoothing of the data, though once again the smoothing is not directly of the 
raw data. If we use a very small T, we are over-smoothing the signal, and we may expect an unreliable 
reconstruction if the actual signal is not smooth enough on timescales < A T/L, where AT is the total 
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length of the data. If on the other hand we use a large L, there are fewer constraints on the Vk, and we 
may get large oscillations between nearby as the fitting procedure tries to compensate for individual 
measurement errors. Using a range of L thus allows an important check on robustness, and in some sense 
checks for systematic errors over and above the statistical errors on the results for a given L. 

When we use a method based on x 2 minimization for finding time delays among a number of curves, 
we may worry about being biased towards large time delays. This is because large time delays imply little 
overlap among the light curves, and so smaller y 2 since there are fewer constraints coming from requiring 
smoothness when the light curves are combined into a single curve. In our method, with a fixed L, this 
bias is partially offset since larger time delays imply larger effective smoothing as the total AT for the 
combined light curves increases. In practice we find negligible bias in the Monte Carlo tests described in 
S3 below. 


We wish to also account for the effect of microlensing in the fitting. Microlensing is the additional 
magnification of each image due to lensing by stars or dark matter clumps in the lensing galaxy. In 
2237+0305 microlensing on the scale of months was first detected by Irwin et al. (1989) and confirmed 
by Houde & Racine (1994). Such microlensing has also been detected in 0957+561, with indications of 
possible microlensing on smaller time scales ( Schild 1996 ). Microlensing may cause an amplification of 
order 0.6 magnitudes, which in PG1115+080 is expected to vary over a time scale ( |Witt et al. 1995| ) of 
~ 10 years (If there is microlensing on a scale of days or weeks in PG1115+080, it appears to be small 
in the data, and we include it in the noise and not explicitly in the analysis). However, even if there is 
no significant variation over the span of the measured light curves, there may be a ’static’ microlensing 
effect which causes different images to vary by different factors in flux, as suggested by SCH. This results 
from the fact that when we observe a quasar varying over time we may effectively be observing two 
different sources. If the flux variation results from a uniform brightening and dimming of the entire 
quasar emission source, then we are only observing a single source. However, if some portions of the 
emission region vary differently from other parts, then the region which produces the mean (i.e. time- 
independent part of the) flux is not identical to the region responsible for the small intensity variations 
with time. A possible physical picture is that the variation is caused by a small jet emanating from the 
core or a sudden localized surge of mass accretion. 


In such a situation, the overall magnifications M t of the mean flux F and m; of the variation / may be 
unequal for a given image i, and their ratio mj/M, may be different for different images. In magnitudes, 
we would then measure (except for an overall factor of —5/2) 


logio (MiF + rruf) 


l°gio (MiF) + log 10 ^1 + -jy-p j 


~ log 10 (MiF) + 


m f 

Mj F log 10 


(3) 


where we assume f/F 1. The variation in magnitudes may therefore have different amplitudes in 
different images if the mi/Mi ratio is not the same for all i. 

The spatial scale for microlensing is set by the Einstein radius of a typical deflector at the lens plane, 
projected onto the source plane. In PG1115+080 this radius is 


Ce = 5 X 10 -3 \(M) /O.IMq] 1/2 pc = 1000 [(M) /O.IMq] 1/2 AU , 


(4) 
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where (M) is the average stellar mass ( |Witt et al. 1995 ). In order to have mj and Mi not be equal, 
one or both of the emission regions must have structure on scales of order (/e or smaller: if they were 
both smooth then the microlensing would be washed out, and we would have r/ij = Mj for all i. On the 
other hand, if the two emission regions overlapped and their union had a very small extent compared to 
(E, then the magnification due to microlensing would not vary over this scale, and again we would have 
m.j = Mj. Thus we reach the interesting conclusion that the indication of different mi/Mi for different 
i (see §3 below) implies that the size of the quasar optical emission region must be of order (e, say to 
within a factor of 10. 


There are, however, a number of reasons to be cautious about this statement. The two emission regions 
may have different sizes, and in this case we might expect the mean region to contain the variation region, 
and then the different mi/Mi would imply the lower bound only on the larger emission region. More 
generally, while there appears to be significant emission from a structure of size C,e , this does not exclude 
a range of structure also on smaller or larger length scales. In addition, </e depends on what we assume 
for (M), and there could be a large range of masses contributing. The intensity variation may involve 
relativistic motions, and such high velocities can cause microlensing variations on small time scales that 
we might confuse as differences in mi/Mi. The stochastic nature of optically-thick microlensing (e.g. 
Witt et al. 1995) implies the possibility of a rare statistical fluctuation over the relatively short time 
span of the current data. Note also that the magnification due to the macrolens itself varies with 
position, but it should only vary over much larger spatial scales. There could be a different amount of 
contamination by light from the lens galaxy in each image, and this would lead to apparent changes in 
the Mj with no corresponding changes in the mj. However, in PG1115+080 this effect seems to be too 
small to affect the mj/Mj ratios significantly. Finally, there may still be unmodelled systematic errors 
in the measured light curves which could give the appearance of unequal mj/Mj ratios. With a longer 
time span of measurements it may be possible to eliminate some or all of these uncertainties. See also 
Gould Sz Miralda-Escude (1997) for an independent, illuminating discussion of the possible observational 
consequences of differential microlensing of accretion disks in quasars. We note that the size given 
in equation ([|) is roughly consistent with the upper limit from the interpretation of microlensing in 
2237+0305 (Wambsganss, Paczynski & Schneider 1990). It is somewhat smaller than an estimate by 
Schild (1996) based on variations in 0957+561. 


In measuring the quasar image light curves, SCH used two nearby stars, *B and *C, as photometric 
references. All intensities (i.e. for the images as well as for *B) were measured with respect to *C, and 
SCH then subtracted one half of the magnitude of *B from the quasar components, for each observation. 
The correlation between the quasar and *B may result from systematic errors such as the adopted point 
spread function, which varies across the held and so may not be identical in *B and in the quasar 
components. In our analysis we treat the fluctuations of *B as stochastic noise, with a component which 
is common to *B and to all the quasar components, as well as an independent source of error. We then 
let the x 2 minimization decide the size of the error common to *B and all the quasar images, separately 
in each observation, but guided by our best estimate of the typical size of the correlation. Thus, we add 
*B to our method effectively as a fourth light curve, but for the underlying model of *B we take just a 
constant intensity in time. We also assume an error component common just to the three quasar images. 
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Note that such correlated errors are visually apparent from the data, e.g. in Figure 2 of SCH near days 
100, 130, and 150, with the clearest example near day 90. 


3. Results for PG1115+080 


Our model for component A consists of L points as described in §2. For B we use the same model, but 
with an overall magnitude shift ms a of A with respect to B (accounting for mean magnification), a time 
shift tsA due to the time delay, and an overall variation ratio ckbA multiplying the variation (accounting 
for differential microlensing as explained in §2). Similarly, we have tuaCi t-AC , and olac-, as well as msc, 
tsci an d olbc ■ When we fit all the curves simultaneously, we have the constraints rriBC = bib A + biac, 
tsc = tsA + tAC, and olbc = olbaolac- The only parameters which require non-linear % 2 minimization 
are tsA, tAC > olba, and olac • We minimize with respect to these parameters using direction set methods 
(press et al. 1992c, §10.5). 


We vary a number of parameters in the model. We try L = 10, 20 and 30 interpolation points. We try 
linear, cubic and cubic spline interpolation. For the measurement errors, we begin with the uncorrelated 
errors for A, B and C of 2.3, 7.5 and 4.0 millimagnitudes (mmag) at lcr, respectively (SCH). We also take 
an uncorrelated error of 1.2 mmag for the bright star *B. We then add an error of standard deviation 
<r n common to *B and all QSO components, and also separate errors of f±a n in *B and f 2 &n hr all the 
QSO components. If we assume /i = /2 = 1, then since *B has a measured variability of 8 mmag at lcr, 
we take cr n — 6 mmag. As we show below, reasonable variations in o n , f\ and fi affect the minimum 
X 2 value but not appreciably values of the derived parameters such as the time delays. On the other 
hand, having no correlated errors at all (a n = 0) does change the derived parameters. All the errors are 
naturally assumed to be Gaussian in the x 2 method. 


We can also take Gaussian errors for Monte Carlo trials. On the other hand, since the method 
reconstructs a model for the input signal, it also reconstructs an estimate of the measurement errors in 
all the light curves in each observation. We take these actual errors instead of Gaussian errors for a 
second, bootstrap-like Monte Carlo analysis. Thus we have a set of 36 observations, and in the Monte 
Carlo trials for each day we pick one of the 36 at random (with replacement) and add the reconstructed 
errors on that day (in A, B, C, and *B) to the matching simulated light curves (i.e. A, B, C, and *B, 
respectively). Note that this is not a rigorous bootstrap, since after the x 2 minimization the errors are 
no longer independent. However, this procedure should give us an idea of the actual size and mutual 
correlations of the measurement errors, and should be more robust and conservative than the Gaussian 
Monte Carlo trials. For the input signal in the Monte Carlo trials we take a smooth curve (see Figure 1) 
with roughly the same shape as the variation in the actual data. This way, the input is identical to the 
real signal on large time scales, and we assume that any apparent small scale variability is dominated by 
measurement noise, which seems likely. 


We define our ’standard’ model as having L = 20, cubic spline interpolation, a n 


6 mmag and 





fi = /2 = 1. Thus the assumed covariance matrix of errors for this model in mmag squared is 


77.3 

72 

72 

36 

72 

128.3 

72 

36 

72 

72 

88.3 

36 

36 

36 

36 

73.3 


( 5 ) 


with rows from the top in the order A, B, C, and *B. Hereafter we use these standard settings except 
where otherwise noted. 

We begin by fitting each quasar component separately, always along with *B. Table 1 shows the 
resulting values of the y 2 per degree of freedom y 2 , for various L. The number of data points is 72 (36 
each for a quasar component and for *B), and the number of parameters is L + 1 (1 for the mean value 
of *B), which yields 71 — L for the number of degrees of freedom (ndof). The reduced y 2 is in most cases 
less than 1, which is reasonable since we are assuming that a significant portion of the error is correlated 
among all the quasar components, and if we take each of them separately then we are not including this 
very strong constraint. Table 1 also suggests that if our error estimates are roughly correct for A and 
C, they may be somewhat underestimated for component B. Note that the y 2 value is higher than 1 at 
L = 10 only for B, which is the faintest component and so has the largest errors. This may indicate 
higher non-Gaussianity for the errors in this component, since with lower L the fewer parameters cannot 
effectively compensate for outlier points. 

Next we fit *B and two quasar components at a time, with the results displayed in Table 2. There 
are 36 * 3 = 108 data points, and L + 4 parameters (e.g. *B mean, msc , t bc and asc)- Also shown are 
the one sided 68% confidence limits derived from 300 bootstrap Monte Carlo trials. Both the bootstrap 
errors and the parameters (tsc etc.) that we input into the Monte Carlo trials are in each case (i.e. for 
each pair of components and for each L) set according to the values fitted from the data in that case. 
So e.g., for the BC pair with L = 20, a value of tsc = 25.3 days is used as input for the Monte Carlo 
trials which in turn determine the error of 1 2 q days on tsc- We only show the results of varying L since 
this leads to the largest variations (see below). Note that the BA time delay has the largest fractional 
uncertainty, as we expect: B has the largest errors, so the uncertainty in tsA is greater than in tAC, 
and while C has larger errors than A, the BC delay is larger and much better resolved than is tsA- For 
the variation ratios, there is some variation with L, but only q.bc is consistent with 1, while the others 
disagree with 1 at the 4 — 5cr level for each L. 

Finally we combine all three quasar components, with results shown in Table 3. Uncertainties, where 
shown, are again one sided 68% confidence limits derived from bootstrap Monte Carlo trials. All the 
models are based on our standard inputs, except for the changes shown in the first column. In rows 4 
and 5, a n is in mmag. In rows 6 and 7, f\ and fa refer to assumptions about the correlated errors (see 
the beginning of this section). In rows 8 and 9, we use cubic and linear interpolation, respectively (see 
§2 above). As above, tabc = tAC /t BA- We try a large range of model assumptions, most of which show 
very little variation in parameter values relative to varying L. The exception is the last row, marked 
’SCH-like’. Here we use our method, but with assumptions that correspond closely to those of SCH: we 
assume no correlated errors, the uncorrelated errors are doubled, *B is not included in the fitting but 
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rather we subtract half of *B from the quasar light curves, and we set the variation ratios to 1. The 
’SCH-like’ parameter values are compatible with those of SCH to about lcr. 

A y 2 method, in which errors are assumed to be Gaussian, can be sensitive to a small number of 
outliers. We use two procedures to check for the effects of such ’bad’ measurements. First, we iteratively 
remove outliers. After fitting with our standard model parameters, we compare each point on every light 
curve to the reconstructed signal. We divide the difference by the individual uncorrelated error for that 
point, which yields the number of standard deviations that each point is off by. We add 20% to the error 
for points off by more than 4a, and add 10% for points off by 2 a — 4a. We then fit with the adjusted 
errors, and repeat until convergence. The resulting values are tsc = 24.0 days and rABC = 1-15, well 
within la of the results from the standard fit. As a second check, we try fitting the time delays after 
removing points from the data set. We remove three points at a time, i.e. one point from each light curve, 
all taken in the same observation. This yields 36 sets of data. In two of these data sets we remove points 
from the first or last observation, thereby changing the total time span of the data and thus the model 
spacing for a fixed L. In order to avoid this, and to be able to check the effects of outliers independently 
from the effects of model spacing, we in fact do not explicitly remove any points, but rather multiply their 
individual errors by 1000 for a de facto removal. The 36 fits to these data sets yield a mean tsc — 25.1 
days with a standard deviation of 0.35 days, a minimum of 23.9 days and a maximum of 26.1 days. They 
also yield a mean rABC = 1.13 with a standard deviation of 0.03, a minimum of 1.05 and a maximum of 
1.20. Our two checks suggest that the fitting results are not strongly influenced by outlier points in the 
data. 

For the standard case ( L = 20), there are 4 * 36 = 144 data points, and L + 7 parameters. In each 
case, we can calculate czbc = ct-BA&AC- E.g. for L = 20 the result is asc = l-lOl^i- With L = 20, the 
lcr uncertainties from Gaussian Monte Carlo trials are, for comparison with Table 3, ]to'g days for tAC > 
days for t B A, ±0 days for t BC, t'H for r A BC, for a AC, i ;°4 for &BA, and for a B c- The 
95% confidence limits from bootstrap Monte Carlo trials are, also for the L = 20 case, ^3 2 days for Cics 
+I 2 days for t B A , days for tsc, tii for r ABC 1 tio for a AC, ^12 for a BA, and 1;^ for a B c■ For this 
standard case, Figure 1 shows the quasar components A (*), B (•), and C (o) (including the final point 
in each light curve, which was excluded from the fitting). All have been corrected in each observation 
only by the error component which is common to all the quasar components and to *B, as reconstructed 
by the fitting. Other than this A is shown as observed (except for a vertical shift to have zero mean), 
but B and C are shifted and scaled according to the values of m.AC-, bib A, tAC-, tsA, otAC: and as a- Also 
shown are the reconstructed signal (solid line) and the input used in the Monte Carlo trials (dotted line). 
Figures 2 and 3 show one dimensional y 2 plots, as a function of tsc and vabc , respectively. In each plot, 
at every point the parameter shown is fixed and y 2 is minimized with respect to all other parameters. 
The plots do not show any strong local minima that could be confused with the global minimum. The 
formal la uncertainties derived from these curves are 1.0 day for tsc and .12 for tabc■ Figure 4 shows 
a two dimensional y 2 plot as a function of tAC and tsA > around the minimum (marked x). The same 
y 2 contours are also shown as a function of tsc and rABC- Note that the result of SCH (marked o) 
is outside even the outermost contour, which delineates the formal 99.99% confidence level. The figure 
clearly shows that the main disagreement of the present results with SCH is in the value of rABC-, while 
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the values of tsc are consistent with each other. 


We adopt the L = 20 results and bootstrap uncertainties, since they are compatible with the values 
for the other L to about la. Thus tsc = 25.0lJ'| days, and Tables 2 and 3 indicate that this 6% 
uncertainty (at 68% confidence) is a reasonable estimate. Our value used with lens models reduces the 
induced Hq only by 5% relative to the SCH value. On the other hand, tabc varies fairly strongly with 
different assumptions, and is only weakly constrained at tabc = 1 • 13lt if • This result is close to the 
values around 1.4 predicted by lens models, but the uncertainty is too large to be able to decide among 
different types of models ( Keeton &; Kochanek 199^ ). The fitting also recovers biac = 2.033 ± .005, 
tuba = —2.534 ± .007, and msc = —-501 ± .008, which yield variation-subtracted magnification ratios 
for the mean flux. These agree with SCH and also with the flux ratios of Kristian et al. (1993), since the 
quasar does not vary much over the time scale of the time delays. Our extremely accurate magnification 
ratios are not useful for lens modelling, since they are still likely to be greatly altered by microlensing, 
as suggested by the variation ratios. 


Regarding the variation ratios, olbc = 1-lO^Qg is consistent with 1, while o.ac = l-39%gg is greater 
than 1 and asA = -79^ Qg is less than 1, both at 4 — 5cr. These values are also roughly consistent with 
Table 2. To make perhaps a more direct test of the significance of this result, we also perform Monte 
Carlo tests on input which has variation ratios equal to 1, but they are allowed to vary in the fitting. The 
result indicates that if the variation ratios were really 1, we would measure olac = l - OOl'gy iX'yj a t 95 % 
confidence) and asA = 1-00(1 23 at 95% confidence), so in this sense the L = 20 results for chac and 
as a are significantly different from 1 at 5cr and 2cr, respectively. As discussed in §2, this may indicate 
that the quasar optical emission region must be of order 1000 AU in linear dimension. 


4. Conclusions 

We have developed a method for determining time delays among light curves of multiple images 
of a gravitational lens. The method constructs a simple model for the actual source variation, using 
interpolation between a number of equally spaced values. It then performs a combined y 2 minimization 
by fitting all of the light curves to this model simultaneously, which is similar to the method of Press et al. 
(1992a, b). The ability to vary the basic parameters of the model over a large range lends robustness to 
our method. Most of the parameters are linear and so the y; 2 minimization is easily done. In addition to 
the time delays, the other non-linear parameters are relative variation ratios, which account for different 
fractional variation in flux for different images. We interpret this physically as evidence for differential 
microlensing, i.e. a different magnification due to microlensing for the varying region from the region 
giving rise to the mean flux. 

Applying our method to the light curves of PG1115+080 observed by SCH, we find a value of 25.0^3^ 
days (95% confidence) for the delay between components B and C, and a ratio tAc/^BA for the two 
smaller delays of 1.131'}* (68% confidence). Unlike SCH, we include correlated measurement errors as 
well as the above mentioned variation ratios in the analysis. Our result for tsc agrees with SCH, but the 
ratio tabc does not. Our result for tabc does agree with lens models, but we find that with the present 
data it cannot be derived accurately enough to help in fitting lens models. For the variation ratios, we 
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find aAC = 1.39l;2o an d «BA = -79^'(95% confidence), each indicating differential microlensing at a 
significance of 4 — 5 times our estimated lcr uncertainties. If confirmed as the data accumulates, this 
would imply that the size of the quasar optical emission region is of order 1000 AU, for microlenses of 
(M) = O.IMq. Further data may also allow a determination of the time variation of the two microlensing 
magnifications. 

I thank Paul Schechter for suggesting this project and for helpful advice, Ed Bertschinger for encour¬ 
agement and useful comments, and Doug Richstone and Ned Wright for valuable comments. This work 
was supported by NASA grant NAG5-2816 and NSF grant AST-9529154. 
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Table 1. Single quasar component fitting 


L 

ndof 

A, X 2 

B, X 2 

c, x 2 

10 

61 

0.66 

1.15 

0.66 

20 

51 

0.58 

0.78 

0.58 

30 

41 

0.67 

0.74 

0.67 



Table 2. 

Fitting of quasar components in pairs 


L 

ndof 

X 2 time delay [days] 

var ratio 

BC components 

10 

94 

1.77 

25-0±l -l 


20 

84 

1.54 

oc q+2.3 
zo.o _2 o 

i-09l:88 

30 

74 

1.03 

26.7iJ;| 

i.03i:88 

AC components 

10 

94 

1.27 

11 n + 0 - 9 

i S5+- 06 

20 

84 

1.06 

n-o+8:l 

i CO +-09 
i - oz -.08 

30 

74 

0.98 

12-3±i;g 

1.401:88 

BA components 

10 

94 

2.38 


o.74i:88 

20 

84 

1.48 

6-751^ 

0.61181 

30 

74 

1.06 

7.22lJ;| 

o.46i:88 
























Table 3. 


Fitting of the three quasar components 


Input* 

ndof 

x 2 

t A c [days] 

tBA [days] 

t-BC [days] 

r ABC 

CVAC 

a BA 


L = 10 

127 

1.83 

1 9 ^+0.9 


26.5111 

0.891:11 

1 4^+- 06 
i - 4O -.06 

,8ol:l 

L = 20 

117 

1.60 

i3.3±?;g 

n-7±l:| 

25.0111 

1.13l;$ 

1 Q O+-07 
i ' oy -.08 

. 791:1 

L = 30 

107 

1.38 

I4.9il;f 

io.6+l:| 

25.5111 

1411:1 

1 4Q+-08 

79+* 04 

* ‘ .05 

&n = 4 

117 

2.27 

13.0 

11.8 

24.8 

1.10 

1.38 

.80 

a n = 8 

117 

1.34 

13.4 

11.7 

25.1 

1.15 

1.39 

.79 

/ i =0 

117 

2.01 

12.8 

11.4 

24.2 

1.12 

1.34 

.80 

/ 2 = 0 

117 

2.13 

12.8 

11.8 

24.6 

1.08 

1.38 

.81 

cubic 

117 

1.56 

13.4 

11.9 

25.2 

1.13 

1.39 

.80 

linear 

117 

1.55 

13.5 

11.7 

25.1 

1.15 

1.40 

.79 

SCH-like 

84 

1.61 

8-9 t\l 

ii-9±l;l 

20.9l|l 

0 . 751:1 

1 

1 


See text. 
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time [ days ] 

Fig. 1.— Light curves for components A (*), B (•), and C (o). B and C have been shifted to match A, 
and all have been partly corrected for errors (see text). Also shown are the reconstructed signal (solid 
line), the input used in Monte Carlo trials (dotted line), and the final point in each light curve, which 
was excluded from the fitting. 



Fig. 2.— Value of x 2 versus assumed delay tsCi relative to the x 2 value at its global minimum. 
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Fig. 4.— x 2 contours in the tAC , tsA plane (left) and in the tsc, 1'ABC plane (right). The point marked 
x is x 2 = 187, the global mimimum. The contours are drawn at Ay 2 =2.30, 4.61, 6.17, 9.21, 11.8 and 
18.4, the la, 90%, 2cr, 99%, 3er and 99.99% confidence levels for two parameters. The point marked o is 
the SCH result. 






